Optimising CNT-FET biosensor design through modelling of biomolecular electrostatic gating and its application to β-lactamase detection

Carbon nanotube field effect transistors (CNT-FET) hold great promise as next generation miniaturised biosensors. One bottleneck is modelling how proteins, with their distinctive electrostatic surfaces, interact with the CNT-FET to modulate conductance. Using advanced sampling molecular dynamics combined with non-canonical amino acid chemistry, we model protein electrostatic potential imparted on single walled CNTs (SWCNTs). We focus on using β-lactamase binding protein (BLIP2) as the receptor as it binds the antibiotic degrading enzymes, β-lactamases (BLs). BLIP2 is attached via the single selected residue to SWCNTs using genetically encoded phenyl azide photochemistry. Our devices detect two different BLs, TEM-1 and KPC-2, with each BL generating distinct conductance profiles due to their differing surface electrostatic profiles. Changes in conductance match the model electrostatic profile sampled by the SWCNTs on BL binding. Thus, our modelling approach combined with residue-specific receptor attachment could provide a general approach for systematic CNT-FET biosensor construction.


Attachment parameterization for H-REMD simulations
The covalent attachment between BLIP2 and the nanotube is represented with an atomistic resolution to precisely consider the constraints imposed by the attachment site on BLIP2's orientations with respect to the nanotube during the Hamiltonian replica-exchange molecular dynamics (H-REMD) simulations.
The BLIP2 azF variants were covalently attached to the single walled carbon nanotube (SWCNT) in the biosensor using azido-phenylalanine photochemistry 5 .In our simulations, the attachment is modelled by a 4-azido-L-phenylalanine, in which two of the three nitrogen atoms of the azido group were removed (Supplementary Figure 11).The remaining nitrogen was placed 0.126 nm from the nanotube and the plane of the phenyl group was oriented perpendicularly to the nanotube surface, as is expected from the known [2+1] cycloaddition reaction, whereby the nitrene radical formed on exposure to UV light inserts perpendicular to the SWCNT 6 .Finally, the perpendicular attachment with the SWCNT was maintained during the simulations by freezing the nitrogen atom and the last carbon of the phenyl ring, while all the other atoms of the attachment residue were free to move.
We parametrized the model compound shown in Supplementary Figure 11c.A phenylalanine was capped at its N-terminal using an acetyl group and at its Cterminal using a N-methyl amide group, as done when the partial charges of amino acids were parametrized for AMBER.Moreover, a N(CH3)2 group was added to the CZ carbon to mimic the attachment with the nanotube, without having to consider the whole carbon nanotube, which is computationally demanding when performing the required ab initio calculations to determine the partial charges.
The bonded and Lennard-Jones parameters of the attachment were determined from the Generalized AMBER forcefield (GAFF) 7 using AmberTools20 8 and ACPYPE 4 to produce GROMACS compatible files.The atom types determined for the model compound are displayed in Supplementary Figure 11d.As the remaining nitrogen atom of the phenyl azide (atom 19) and the last carbon of the phenyl ring (atom 17) are frozen during the simulation, we do not explicitly define bonded interactions between the nanotube and the attachment.
The partial charges of the atoms were determined using the restrained electrostatic potential (RESP) protocol as done to optimize the partial charges of AMBER's amino acids 9 .We considered two backbone structures for the model compound: one representative of a helical state (=−60° =−60°, =−40° =−40° and =180° =180°) and the other of an extended state (=217° =217°, =160° =160° and =180° =180°), as done for phenylalanine in AMBER.First, we determined the electrostatic potential generated by the two conformations using Gaussian16 10 to perform ab initio Hartree-Fock calculations at the 6-31G* level, as done for AMBER.Second, the RESP protocol was applied simultaneously on the two conformations to optimize the partial charge to best represent the electrostatic potential determined using ab initio calculations.RESP follows a two-step procedure: (1) the charges of all atoms are optimized, except those that are constrained, then (2) the charges of the CH3 groups are reoptimized.During RESP, we constrained the partial charges of the N-terminal acetyl group and the C-terminal N-methyl amide group as well as the partial charges of the mainchain C=O and N-H groups to those in AMBER.Moreover, the total charge of the model compound is constrained to zero and chemically equivalent atoms are constrained to the same charge.The optimized charges are displayed in Supplementary Figure 11d.The electrostatic potential generated by the optimized charges has a relative root mean square of 0.12 with respect to the ab initio electrostatic potential, which is similar to the values obtained for AMBER 9 .

H-REMD simulations
We performed Hamiltonian replica-exchange molecular dynamics (H-REMD) simulations to exhaustively sample the orientations of BLIP2 with respect to the carbon nanotube, while considering the interactions between them.The goal of H-REMD simulations is to enhance the thermodynamics sampling of a molecular system by decreasing the energetic barriers that impair the sampling of traditional MD simulations 11 .To do so, several MD simulations are simultaneously launched in parallel and ordered as followed: the Hamiltonian (energy) in the first MD is unchanged, while it is progressively changed in the other MDs to favour transitions over the energetic barriers.Exchanges of the system between neighbouring MDs is periodically tried using the Metropolis criteria to allow energetically favourable states to diffuse back towards the first MD, which is described by the physically relevant, unscaled energy.
For our system, we realized that BLIP2 hardly samples different orientations with respect to the nanotube when using standard MD simulations.After a quick initial reorganization of a few nanoseconds, the orientation of BLIP2 remains mostly constant for the remaining of the 1000-ns MD simulation (Supplementary Figure 12).Moreover, we realized that the temperature replica-exchange molecular dynamics (T-REMD) technique usually used to get exhaustive sampling when characterizing biomolecule-SWCNT interactions isn't appropriate for our system because the required high temperatures to unbind the protein from the nanotube aren't preserving the overall structure of the protein.Therefore, we choose to use the H-REMD technique in which we selectively scale the attractive part of the Lennard-Jones interactions between BLIP2 and the carbon nanotube.As a result, BLIP2 dissociates more easily from the nanotube at larger scaling, allowing it to thoroughly sample various orientations with respect to the nanotube as the system diffuses back to lower scaling, while preserving the overall structure of BLIP2.
One H-REMD simulation is performed in explicit solvent for each attachment site.The H-REMD simulations were performed using the GROMACS software 2,12 version 2019.6 augmented with the open-source, community-developed PLUMED library version 2.6.2 13,14 for H-REMD simulations 15 .Each H-REMD simulation consists of 24 MD simulations in parallel with scales ranging from 1.0 to 0 that are applied on the Lennard-Jones attractive interactions between BLIP2 and the nanotube.The systems are solvated and equilibrated as follow: (1) the systems are solvated in explicit water molecules and in 150 mM NaCl to mimic the 1xPBS experimental condition, (2) an additional eleven Na ions are added to neutralize the systems, (3) the systems are minimized first using steepest descent and then conjugate gradient, (4) the temperature is equilibrated at 300 K in a 1-ns simulation during which the position of the heavy atoms of the protein is restrained, and (5) the pressure is equilibrated semi-isotropically, except on the direction along the nanotube axis, at 1 atm using the Berendsen barostat in a 1-ns simulation during which the position of the heavy atoms of the protein is again restrained.
The initial systems for the H-REMD simulations are shown in (Supplementary Figure 13).To generate the initial states of the scales, we performed a 50-ns MD simulation during which the attractive part of the Lennard-Jones interactions between BLIP2 and the nanotube were removed (scaling of 0) to generate various orientations of BLIP2 with respect to the nanotube, without any contact between BLIP2 and the nanotube.From that simulation, we clustered the trajectory using a cutoff of 1 nm on the carbon alpha atoms of BLIP2 (without aligning the structures since the nanotube is fixed).We used the clusters' centers as initial configurations for the scales to ensure a representativity of the various orientations accessible to BLIP2 when it does not interact with the nanotube.

Electrostatic potential simulations
The APBS software was used to solve the non-linear Poisson-Boltzmann equation 16 to determine the electrostatic potential generated by the protein on the SWCNT.In these simulations, the van der Waals radius and the partial charge of each atom were set to those of AMBER14sb, while the partial charges of the linker atoms were set to 0 to focus on the electrostatic potential generated by the protein.The relative permittivity of the solvent and non-solvent were respectively set to 78.54 and 2. The solvent region was determined using a water molecule radius of 0.140 nm to determine the solvent boundary of the BLIP2/attachment/nanotube system.A monovalent salt concentration of 150 mM was used to mimic the PBS experimental conditions and the radius of the ions was set to 0.2 nm.The electrostatic potential values on the boundary of the simulation box were determined using the simple Debye-Hückel model.
The electrostatic potential was determined using a two-step procedure: the nonlinear Poisson-Boltzmann equation was first solved on a large box having 1000 points/nm 3 (1 point/Å 3 ) so as to consider the system in its entirety with a large enough space system and boundaries, then the non-linear Poisson-Boltzmann equation was again solved to refine the calculated electrostatic potential values on a region of 2 nm around the nanotube using 8000 points/nm 3 (8 points/Å 3 ).In the first step, the size of simulation box was determined from the size of the BLIP2/TEM-1 complex in the simulations multiplied by 1.7 along each direction, yielding a box of 22.5 nm by 22.5 nm by 16.1 nm.The center of the box was set to the average center of the protein in all the sampled configurations.In the second step, the region of 2 nm around the nanotube corresponds to a box of 17.6 nm by 4.80 nm by 4.80 nm that was centred on the nanotube.The length of the nanotube was doubled compared to the H-REMD simulations to consider the fixed Debye-Hückel boundary conditions used in the electrostatic potential calculations.This two-step procedure was necessary to accurately solve, in a reasonable time the electrostatic potential of the thousands of configurations generated by the H-REMD simulations.The two-step procedure yields similar electrostatic potential values on the nanotube compared to only doing the first step using the greater resolution of the second step (8000 points/nm 3 ).The difference between the two approaches is only 0.01 mV on average, which is small compared to the magnitude of the electrostatic potential on the carbon nanotube (<100 mV).
We considered all configurations sampled at all Hamiltonian scales for the analysis on the 500 to 1250 ns interval with a 10-ps resolution (75000 configurations per scale).To do so, we determined the normalized weight at the unscaled energy of each configuration using the weighted histogram analysis method (WHAM) 17 .For the electrostatic potential, we use the configurations sampled over the 500 to 1250 ns interval with a 100-ps resolution (7500 configurations per scale).Also, we focused the analysis on the configurations sampled in the Hamiltonian scales 1 to 7 because they have the greatest normalized weights (>10 -9 compared to a maximum of 10 -5 ).The AF models can be found here: https://github.com/drdafyddj/CNT-electrostaticmodelling.
Supplementary Figure 3. Rotamers azF that form steric clashes with an SWCNT on manual docking.The red and blue residues represent rotamer 1 of BLIP2 41azF and rotamer 3 of BLIP2 213azF , respectively, with yellow residues sterically clashing with the nanotube.Docking of BLIP2 azF mutant onto the SWCNT is based on the known nitrene insertion chemistry reaction to provide initial models of the CNT-receptor protein configuration.The nitrene radical formed on exposure to UV light is expected to insert perpendicular to the SWCNT 6 .Attachment of BLIP2 41azF or BLIP2 213azF resulted in an averaged height increase of 1.5 nm or 1.9 nm, respectively.This is lower than the heights anticipated from the dimensions of BLIP2 but is in line with previous AFM measurements that routinely underestimate the heights of soft materials such as proteins [23][24][25] .After exposure to TEM-1, the height of the CNT-FET increased slightly from 1.5 nm to 2.1 nm and 1.9 nm to 2.3 nm for BLIP2 41azF and BLIP2 213azF , respectively.These relatively small increases (0.4 -0.6 nm) potentially provides an indication on the position of BLIP2's β-lactamase binding face; with adjacent binding to the SWCNT more likely to fit this data than the β-lactamase binding above the plane of the SWCNT.Visualising single molecule attachment AFM provides a valuable insight to peripheral SWCNTs that branch off larger bundles.These small-scale nanotubes offer potential to visualise single molecule attachment, with the ~2-5 nm sized proteins notable against single SWCNTs.Supplementary Figure 3d above shows the clearest example of this, with clear height changes (see white arrows) at the end of a SWCNT after BLIP2 41azF attachment; the round shape is consistent with that of a protein.Considering the AFM image of the NT-FET was taken after a washing step following protein attachment (see methods), it can assumed that these species are BLIP2 41azF rather than adsorbed contaminants.The same SWCNT was then located after the sensing experiment with KPC-2, where initial protein areas appear to have expanded.This suggests KPC-2 has bound adjacently to BLIP2 41azF .The extracted height profiles (Supplementary Figure 3e) support this hypothesis, with the defined protein peaks from BLIP2 41azF functionalisation broadening and increasing in height after complexing with KPC-2.This provides further proof that BLIP2 retains BL binding function on SWCNT attachment and provides and interesting insight into protein binding morphology.

Supplementary
scales in the 500 to 1250 ns interval, where scale 1 (purple) is the unscaled energy and scale 24 (red) has no attractive protein-CNT interaction.The accessible orientations were determined from rigid body rotations of the attached protein using three dihedral angles bridging azF and the backbone, namely the dihedrals 14-12-11-8, 12-11-8-7 and 11-8-7-1 illustrated in Supplementary Figure 11a-c.From these orientations, the minimal distance of any atom of the protein to the CNT was determined to identify clashes with the CNT, where 0 nm (black) means an inaccessible orientation unless structural rearrangement in the protein.Source data for plots b and c can be found here: https://zenodo.org/records/12783255.Supplementary Discussion.In panel c, we observe that the H-REMD simulations sample most of the accessible orientations, including most of those that are near the CNT (the white-to-black transition region).As expected, the orientations sampled at low scaling (scales 12 and smaller) are those for which the protein is nearer the CNT, while the orientations sampled at high scaling (scales 12 and greater) are those for which the protein is not in contact with the CNT because the attractive part of the Lennard-Jones interaction between the protein and the CNT are progressively scaled in the H-REMD simulations by a factor going from 1 (unscaled: scale number 1) to 0 (completely scaled: scale number 24).In panel b, we observe the most energetically favorable orientations at the unscaled energy (scale 1), corresponding to a subclass of the orientations for which the protein is more tightly in contact with the CNT.
the azF mutation within the binding interface and thus will block BL binding on attachment to a SWCNT.(c) Concentration-dependent change in conductance of BLIP2 213azF functionalised devices on addition of non-specific protein, Bovine serum albumen (BSA).(c) Concentration-dependent change in conductance of BLIP2 213azF functionalised devices on addition of TEM-1 in the presence of 0.5% (w/v) BSA.The conductance data were collected via I-V measurements between source and drain In all cases data is presented as the mean value with errors bars representing +/-the SD (n=3 where n is the number of different devices measurements were derived from).For S10a-c no fit to the data was possible so trendlines are shown as dashed lines.Source data for plots are provided in the Source Data file.

Figure 4 .
AFM analysis of the CNT-FET transistor channel.AFM images on the left show the functionalised transistor channel, with the overlaid horizontal white lines indicating six cross sections taken for height analysis.An average cross section was calculated for the channel before and after BLIP2 azF functionalisation, and again following β-lactamase sensing, and are plotted in the graphs on the right.Height profiles are seen for SWCNTs (black dashed), BLIP2 azF ), TEM-1 (blue).(a) BLIP2 41azF derivatised CNT-FET, (b) BLIP2 213azF derivatised CNT-FET and (c) relative changes in height of the base SWCNT on attachment of BLIP2 and binding of TEM-1.(d) AFM analysis of single/double SWCNT bundles branching off the main transistor channel reveal individual BLIP2 41azF proteins attaching along the length of it.Following the sensing experiment, KPC-2 proteins appear to bind adjacently.(e) Extracted height profile before and after BLIP2 41azF functionalisation, and after KPC-2 sensing.Arrows identify speculative protein peaks.Source data for plots in panels a, b and e are provided in the Source Data file.Supplementary Discussion.